#########################################################################################
# R script for
# Tober, T. & Busemeyer, Marius R.
# 'Breaking the Link? How European Integration Shapes Social Policy Demand and Supply'
# Journal of European Public Policy
# Table 1 and Figure 2
#########################################################################################

#########################################################################################
# Load packages
#########################################################################################

if (!require("brms")) install.packages("brms")
library(brms)

#########################################################################################
# Load data
#########################################################################################

load("data_mixed.RData")

#########################################################################################
# Table 1
#########################################################################################

Sys.getenv('PATH')
system('g++ -v')
system('where make')

# Set priors
prior <- c(set_prior("normal(0,1)", class = "b"),
           set_prior("student_t(4,0,1)", class = "sd")) 

# Model 1
fit0 <- brm(red_pref_log2 ~ agea + gndr + eduyrs + edctn 
            + pdwrk + uemp + rlgdgr + union_mb + lrscale + hinc_dum 
            + sm_mean + sm_diff + conf_mean + cf_diff  
            + (1 | Country) + (1 | Year) + (1 | Country:Year),
            data = datsc, family = bernoulli(link = "logit"), prior = prior,
            warmup = 500, chains = 2, iter = 2000, core = 2,
            control = list(adapt_delta = 0.95))
summary(fit0)

# Model 2
fit1 <- brm(red_pref_log2 ~ agea + gndr + eduyrs + edctn 
            + pdwrk + uemp + rlgdgr + union_mb + lrscale + hinc_dum 
            + sm_mean + sm_diff + conf_mean + cf_diff + socprot_mean + socprot_diff 
            + (1 | Country) + (1 | Year) + (1 | Country:Year),
            data = datsc, family = bernoulli(link = "logit"), prior = prior,
            warmup = 500, chains = 2, iter = 2000, core = 2,
            control = list(adapt_delta = 0.95))
summary(fit1)

# Model 3
fit2 <- brm(red_pref_log2 ~ agea + gndr + eduyrs + edctn 
            + pdwrk + uemp + rlgdgr + union_mb + lrscale + hinc_dum 
            + sm_mean + sm_diff + conf_mean + cf_diff + gdppc_mean + gdppc_diff 
            + (1 | Country) + (1 | Year) + (1 | Country:Year),
            data = datsc, family = bernoulli(link = "logit"), prior = prior,
            warmup = 500, chains = 2, iter = 2000, core = 2,
            control = list(adapt_delta = 0.99))
summary(fit2)

# Model 4
fit3 <- brm(red_pref_log2 ~ agea + gndr + eduyrs + edctn 
            + pdwrk + uemp + rlgdgr + union_mb + lrscale + hinc_dum 
            + sm_mean + sm_diff + conf_mean + cf_diff + uempl_mean + uempl_diff 
            + (1 | Country) + (1 | Year) + (1 | Country:Year),
            data = datsc, family = bernoulli(link = "logit"), prior = prior,
            warmup = 500, chains = 2, iter = 2000, core = 2,
            control = list(adapt_delta = 0.99))
summary(fit3)

# Model 5
fit0.sub <- brm(red_pref_log2 ~ agea + gndr + eduyrs + edctn 
                + pdwrk + uemp + rlgdgr + union_mb + lrscale + hinc_dum 
                + smopen_mean + smopen_diff + smimp_mean + smimp_diff 
                + confpart_mean + cfpart_diff + confcomp_mean + cfcomp_diff
                + (1 | Country) + (1 | Year) + (1 | Country:Year),
                data = datsc, family = bernoulli(link = "logit"), prior = prior,
                warmup = 500, chains = 2, iter = 2000, core = 2,
                control = list(adapt_delta = 0.99))
summary(fit0.sub)

#########################################################################################
# Figure 2
#########################################################################################

tmpdat <- datsc[, c("agea", "gndr", "chldhm", "eduyrs", "edctn", "pdwrk", "uemp", 
                    "rlgdgr", "union_mb", "lrscale", "hinc_dum", "sm_mean", "sm_diff", 
                    "conf_mean", "cf_diff", "Country", "Year")]
jvalues <- with(datsc, seq(from = min(sm_diff, na.rm = T), to = max(sm_diff, na.rm = T), length.out = 10))
tmpdat$sm_diff <- jvalues[1]
predprob1 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[2]
predprob2 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[3]
predprob3 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[4]
predprob4 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[5]
predprob5 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[6]
predprob6 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[7]
predprob7 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[8]
predprob8 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[9]
predprob9 <- fitted(fit0, newdata = tmpdat)
tmpdat$sm_diff <- jvalues[10]
predprob10 <- fitted(fit0, newdata = tmpdat)
predprob_sm <- rbind(colMeans(predprob1),colMeans(predprob2),colMeans(predprob3),
                  colMeans(predprob4),colMeans(predprob5),colMeans(predprob6),
                  colMeans(predprob7),colMeans(predprob8),colMeans(predprob9),
                  colMeans(predprob10))
plotdat_sm <- as.data.frame(cbind(predprob_sm, jvalues))

tmpdat <- datsc[, c("agea", "gndr", "chldhm", "eduyrs", "edctn", "pdwrk", "uemp", 
                    "rlgdgr", "union_mb", "lrscale", "hinc_dum", "sm_mean", "sm_diff", 
                    "conf_mean", "cf_diff", "Country", "Year")]
jvalues2 <- with(datsc, seq(from = min(cf_diff, na.rm = T), to = max(cf_diff, na.rm = T), length.out = 10))
tmpdat$cf_diff <- jvalues2[1]
predprob1cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[2]
predprob2cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[3]
predprob3cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[4]
predprob4cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[5]
predprob5cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[6]
predprob6cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[7]
predprob7cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[8]
predprob8cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[9]
predprob9cf <- fitted(fit0, newdata = tmpdat)
tmpdat$cf_diff <- jvalues2[10]
predprob10cf <- fitted(fit0, newdata = tmpdat)
predprob_cf <- rbind(colMeans(predprob1cf),colMeans(predprob2cf),colMeans(predprob3cf),
                     colMeans(predprob4cf),colMeans(predprob5cf),colMeans(predprob6cf),
                     colMeans(predprob7cf),colMeans(predprob8cf),colMeans(predprob9cf),
                     colMeans(predprob10cf))
plotdat_cf <- as.data.frame(cbind(predprob_cf, jvalues2))

pdf("marginal.pdf", height = 7, width = 12)
par(mfrow=c(1,2))
plot(plotdat_sm$jvalues, plotdat_sm$Estimate, type = "l", lwd = 2, ylim = c(0.2,0.4),
     ylab = "Probability(Demand for social policy)", 
     xlab = "Economic integration (W)",
     main = "A Within-country economic integration")
grid()
lines(plotdat_sm$jvalues, plotdat_sm$Q2.5, col = "#e34a33", lty = 2, lwd = 2)
lines(plotdat_sm$jvalues, plotdat_sm$Q97.5, col = "#e34a33", lty = 2, lwd = 2)
rug(datsc$sm_diff)
plot(plotdat_cf$jvalues2, plotdat_cf$Estimate, type = "l", lwd = 2, ylim = c(0.2,0.4),
     ylab = "Probability(Demand for social policy)", 
     xlab = "Political integration (W)",
     main = "B Within-country political integration")
grid()
lines(plotdat_cf$jvalues2, plotdat_cf$Q2.5, col = "#e34a33", lty = 2, lwd = 2)
lines(plotdat_cf$jvalues2, plotdat_cf$Q97.5, col = "#e34a33", lty = 2, lwd = 2)
rug(datsc$cf_diff)
dev.off()


